Thalamocortical Connections Mature along a Sensorimotor-Association Axis of Cortical Developmental Heterochronicty

Read in Data for Developmental Characterization

Glasser parcel names

#Glasser regions with corresponding labels, tract names, and mesulam assignments
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name 
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) #create tract variable with format thalamus_$hemi_$region_autotrack, no dashes -

#Glasser regions assigned to mesulam zones
glasser.assignments <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360-mesulam_economo_yeo-assignments.csv")
glasser.assignments <- merge(glasser.assignments, glasser.regions, by = "label", sort = F)
glasser.assignments <- glasser.assignments %>% select(label, tract, orig_parcelname, mesulam_assignment)

S-A axis

S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.assignments, by = c("orig_parcelname"), sort = F)
S.A.axis$SA.axis.bin <- as.numeric(cut2(S.A.axis$SA.axis, g = 5)) #divide the S-A axis into 5 ranked bins, 72 regions per bin

Neurosynth term maps

#Neurosynth z-score maps generated by nimare. Nimare uses multilevel kernel density analysis- Chi-square analysis + FDR-correction to use the same procedure as Neurosynth
neurosynth.terms <- read.csv("/cbica/projects/thalamocortical_development/Maps/neurosynth/atl-glasser360_neurosynth.csv") %>% select(-regionName) %>% select(-timing) #neurosynth term meta-analytic scores for 123 terms present in both NeuroSynth and the Cognitive Atlas, ordered in lh --> rh glasser order
neurosynth.termlist <- names(neurosynth.terms) #list of terms (cue willy wonka line "dont you want to know our names?" "cant imagine why it wouldnt matter")
neurosynth.terms$label[1:180] <- glasser.regions$label[181:360] #lh first
neurosynth.terms$label[181:360] <- glasser.regions$label[1:180] #rh

Spatial permutation (spin) test parcel rotation matrix

perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinning

Dataset-specific diffusion measure spreadsheets

FA.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_FA_finalsample_primary.csv") #generated by sample_construction/PNC/tractmeasures_dfs_PNC.R

FA.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_FA_finalsample_primary.csv") #generated by /sample_construction/HCPD/tractmeasures_dfs_HCPD.R

Dataset-specific tract lists

tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R

tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R

Developmental statistics

setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/") #gam output dir

#list all FA development files in development results directory for accessing. Files were generated by gam_functions/fit_ageGams.R
files <- list.files(getwd(), pattern='FA', ignore.case=T, full.names = F) 

filenames <- c()
#generate variable names to assign files data to 
for(name in files){
  Rname <- gsub('^.{12}|.{4}$', '', name) #remove ".csv" from end of filename and "development_" from start of filename
  filenames <- append(filenames, Rname) #save filenames into a character vector 
}

#read in files and assign to variables
for(i in 1:length(filenames)){
  
  Rfilename <- sprintf("%s",filenames[i]) #save filename as an individual string
  
  if(grepl("gameffects", Rfilename)) {
    x <- read.csv(files[i], header = TRUE) 
    x <- merge(x, glasser.assignments, by = "tract")
    assign(Rfilename, x) 
  }
  
  if(!grepl("temporal", Rfilename) & !grepl("gameffects", Rfilename)) {
    x <- read.csv(files[i], header = TRUE) 
    x <- merge(x, S.A.axis, by = "tract")
    assign(Rfilename, x) 
  }
  
  if(grepl("sexeffects", Rfilename)) {
    x <- read.csv(files[i], header = TRUE) 
    x <- merge(x, glasser.assignments, by = "tract")
    assign(Rfilename, x) 
  }
  
}
rm(x)

A Spectrum of Thalamocortical Developmental Trajectories

Cortex-wide thalamic connectivity developmental profiles

Number of significant developmental effects

#Function to calculate the number and percent of thalamocortical connections showing a significant developmental change in a given measure
sig.effects <- function(measure, atlas, dataset){
  ageeffects.df <- get(sprintf("gameffects_%s_%s_%s", measure, atlas, dataset)) 
  ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 #fdr-corrected significant connections
  sigeffect.totaln <- sum(ageeffects.df$significant) #total number of significant connections 
  sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) #percent of significant connections
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes in %s", dataset, sigeffect.totaln, sigeffect.percent, measure))
}

PNC

sig.effects(measure = "FA", atlas = "glasser", dataset = "pnc")
## [1] "In pnc, 201 thalamocortical connections (90.13 percent) show significant age-related changes in FA"

HCPD

sig.effects(measure = "FA", atlas = "glasser", dataset = "hcpd")
## [1] "In hcpd, 180 thalamocortical connections (78.26 percent) show significant age-related changes in FA"

Number of significant age by sex interactions

#Function to calculate the number and percent of thalamocortical connections showing a significant age*sex interaction
sig.effects.sex <- function(measure, atlas, dataset){
  sexeffects.df <- get(sprintf("sexeffects_%s_%s_%s", measure, atlas, dataset)) 
  sexeffects.df$significant <- p.adjust(sexeffects.df$GAM.int.pvalue, method = c("fdr")) < 0.05 #fdr-corrected 
  sigeffect.totaln <- sum(sexeffects.df$significant) #total number of significant connections 
  sigeffect.percent <- round(sigeffect.totaln/length(sexeffects.df$significant)*100, 2) #percent of significant connections
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age by sex interactions", dataset, sigeffect.totaln, sigeffect.percent))
}

PNC

sig.effects.sex(measure = "FA", atlas = "glasser", dataset = "pnc")
## [1] "In pnc, 0 thalamocortical connections (0 percent) show significant age by sex interactions"

HCPD

sig.effects.sex(measure = "FA", atlas = "glasser", dataset = "hcpd")
## [1] "In hcpd, 4 thalamocortical connections (1.74 percent) show significant age by sex interactions"

Thalamocortical connection GAM smooth functions

PNC

smoothcentered_FA_glasser_pnc.plot <- merge((smoothcentered_FA_glasser_pnc %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_pnc %>% select(tract, GAM.smooth.partialR2)), by = "tract") #zero-centered developmental smooth functions for RH connections + tract-specific partial R2 for plotting

ggplot(smoothcentered_FA_glasser_pnc.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(0, 0.1), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_PNC.pdf", device = "pdf", dpi = 500, width = 2.02, height = 1.65)

HCPD

smoothcentered_FA_glasser_hcpd.plot <- merge((smoothcentered_FA_glasser_hcpd %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_hcpd %>% select(tract, GAM.smooth.partialR2)), by = "tract")

ggplot(smoothcentered_FA_glasser_hcpd.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.02, 0.075), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  scale_y_continuous(limits = c(-0.02, 0.0138)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_HCPD.pdf", device = "pdf", dpi = 500, width = 2.02 , height = 1.65)

Region-specific GAM smooths

#Function for connection-specific GAM spline + derivative plots
connection_trajectory_plot <- function(measure, atlas, dataset, tract_name, display_color, y_ticks){
  
  #Read in participant-level data and model-level fitted values and derivatives 
  participant.data.df <- get(sprintf("%s.%s.%s", measure, atlas, dataset))
  fittedvalues.df <- get(sprintf("fittedvalues_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
  derivatives.df <- get(sprintf("derivatives_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
  derivatives.df$significant.derivative[derivatives.df$significant.derivative == 0] <- NA

  #Connection spline plot with participant-level data
  trajectory.plot <- ggplot(data = participant.data.df, aes(x = age, y = get(tract_name))) +
    geom_point(color = alpha(display_color, 0.5), size = .1) +
    geom_line(data = fittedvalues.df, aes(x = age, y = fitted), color = display_color, linewidth = 1) +
    geom_ribbon(data = fittedvalues.df, aes(x = age, y = fitted,  ymin = lower, ymax = upper), alpha = .8, linetype = 0, fill = display_color) +
    labs(x="\nAge", y=sprintf("%s %s: %s\n", tract_name, measure, toupper(dataset))) +
    theme_classic() +
    scale_x_continuous(breaks = c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
    scale_y_continuous(breaks = y_ticks) + 
    theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 5, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  #Significant derivatives (rate of age-related change) plot 
  derivatives.plot <- ggplot(data = derivatives.df) + 
    geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
    scale_fill_gradient(low = alpha(display_color, 0.2), high = display_color,  na.value = "white") +
    labs(x="\nAge", fill = "Derivative") + 
    theme_classic() +
    scale_x_continuous(breaks=c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
    theme(
        legend.position = "none",
        axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size = 6, family = "Arial", color = c("black"))) 
  
  allplots <- list(trajectory.plot, derivatives.plot)
  tract.plot <- cowplot::plot_grid(rel_heights = c(16, 3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)

  return(tract.plot) 
}

Primary motor cortex (4)

Left 4: S-A axis rank = 16, mesulam assignment = idiotypic

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.4, 0.45, 0.5, 0.55))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               16.06784
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.35, 0.4, 0.45, 0.5))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               15.45184

Lateral parietal (PF)

Left PF: S-A axis rank = 248, mesulam assignment = heteromodal

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1                18.4531
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               16.49456

Lateral prefrontal (8BL)

Left 8BL: S-A axis rank = 342, mesulam assignment = heteromodal

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               21.21106
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               21.91667

Maturational patterns are highly reproducible

Correlation of thalamocortical connection development effect sizes between datasets

Pearson’s R

gameffects.merged <- merge(gameffects_FA_glasser_pnc, gameffects_FA_glasser_hcpd, by = c("tract", "label", "orig_parcelname", "mesulam_assignment"), suffixes = c("_pnc", "_hcpd"))

cor.test(gameffects.merged$GAM.smooth.partialR2_pnc, gameffects.merged$GAM.smooth.partialR2_hcpd)
## 
##  Pearson's product-moment correlation
## 
## data:  gameffects.merged$GAM.smooth.partialR2_pnc and gameffects.merged$GAM.smooth.partialR2_hcpd
## t = 16.04, df = 219, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6677655 0.7902850
## sample estimates:
##       cor 
## 0.7349671

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects.merged, by = c("label", "orig_parcelname", "tract")) #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly  reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$GAM.smooth.partialR2_pnc, y = spin.data$GAM.smooth.partialR2_hcpd, perm.id = perm.id.full, corr.type = "pearson") 
## [1] 0

Correlation plot

ggplot(gameffects.merged, aes(x = GAM.smooth.partialR2_pnc, y = GAM.smooth.partialR2_hcpd)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(-0.01, 0.17)) +
  scale_y_continuous(limits = c(-0.06, 0.15)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.99, height = 1.58)

Correlation of thalamocortical connection age of maturation between datasets

Pearson’s R

cor.test(gameffects.merged$smooth.increase.offset_pnc, gameffects.merged$smooth.increase.offset_hcpd)
## 
##  Pearson's product-moment correlation
## 
## data:  gameffects.merged$smooth.increase.offset_pnc and gameffects.merged$smooth.increase.offset_hcpd
## t = 8.6088, df = 153, p-value = 8.379e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4545288 0.6686745
## sample estimates:
##       cor 
## 0.5712442

Spatial permutation (spin) based p-value

perm.sphere.p(x = spin.data$smooth.increase.offset_pnc, y = spin.data$smooth.increase.offset_hcpd, perm.id = perm.id.full, corr.type = "pearson") 
## [1] 2e-04

Correlation plot

ggplot(gameffects.merged, aes(x = smooth.increase.offset_pnc, y = smooth.increase.offset_hcpd)) +
  geom_point(color = c("#9c3a80"), size = 0.4) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0.05, 0)) +
  scale_y_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0, 1)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Agematuration_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.94, height = 1.58)

Correlation of developmental magnitude and age of maturation within dataset

PNC

Spearman’s rho

cor.test(gameffects_FA_glasser_pnc$GAM.smooth.partialR2, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_pnc$GAM.smooth.partialR2 and gameffects_FA_glasser_pnc$smooth.increase.offset
## S = 234341, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8242399

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0

Correlation plot

ggplot(gameffects_FA_glasser_pnc, aes(x = GAM.smooth.partialR2, y = smooth.increase.offset)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nAge effect (partial R2)", y="Age of maturation (years)\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_Ageofmaturation_PNC.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.58)

HCPD

cor.test(gameffects_FA_glasser_hcpd$GAM.smooth.partialR2, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_hcpd$GAM.smooth.partialR2 and gameffects_FA_glasser_hcpd$smooth.increase.offset
## S = 137678, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8161002

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0

Early and Late Maturing Thalamocortical Connections

Early and late maturing areas

#Identify parcels that fall within first and fourth quartiles of the age of maturation variable
matstage.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
matstage.df$maturation_bracket <- "mid" #set default bracket to mid
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset < quantile(matstage.df$smooth.increase.offset, 0.25, na.rm = T)] <- "early" #set first quartile bracket to early
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset > quantile(matstage.df$smooth.increase.offset, 0.75, na.rm = T)] <- "late" #set fourth quartile bracket to late

matstage.df <- rbind(matstage.df, data.frame(label = "rh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey
matstage.df <- rbind(matstage.df, data.frame(label = "lh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey

ggseg(.data = matstage.df, atlas = "glasser", mapping = aes(fill = maturation_bracket, colour=I("black"), size=I(.06))) + scale_fill_manual(values = c("#FEC22F", alpha("#323280", 0.97), "white", "white" ), na.value = "gray90") + theme(legend.position = "none")
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  side  region label smooth.increase.offset maturation_bracket
##   <chr> <chr> <chr> <chr> <chr>  <chr>                  <dbl> <chr>             
## 1 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L…                   16.9 mid

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Maturationalbrackets_PNC.pdf", device = "pdf", dpi = 500, width = 4.4, height = 4.4)

Neurosynth decoding of age of maturation maps

Compute cognitive term-development map spatial correlations

#Function to calculate the correlation between a neurosynth term z-score map and a thalamocortical development map
developmentmap.neurosynth.decoding <- function(df, dev.metric, term){
  df.allregions <- left_join(spin.parcels, df, by = c("label", "tract", "orig_parcelname")) #full set of parcel data in rh --> lh order 
  term.cor <- cor(df.allregions[,term], df.allregions[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between neurosynth term map and development metric
  term.results <- data.frame("term" = term, "term.correlation" = term.cor)
  return(term.results)
}

Developmental timing decoding by Neurosynth cognitive terms

PNC

gameffects_neurosynth_pnc <- merge(gameffects_FA_glasser_pnc, neurosynth.terms, by = "label")

devpattern.neurosynth.pnc <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_pnc, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
  arrange(desc(term.correlation))
devpattern.neurosynth.pnc <- rbind(slice_max(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map

devpattern.neurosynth.pnc
##                   term term.correlation
## 1           expectancy        0.4099804
## 2           monitoring        0.3968318
## 3    cognitive_control        0.3822346
## 4            reasoning        0.3317204
## 5  response_inhibition        0.3131374
## 6             strategy        0.3119431
## 7     memory_retrieval        0.2883541
## 8              thought        0.2863203
## 9               memory        0.2799527
## 10              recall        0.2692535
## 11  object_recognition       -0.2919385
## 12   visual_perception       -0.2896572
## 13          morphology       -0.2747329
## 14       consolidation       -0.2571732
## 15                gaze       -0.2507995
## 16           expertise       -0.2367694
## 17        coordination       -0.2367162
## 18          perception       -0.2249922
## 19   facial_expression       -0.2183252
## 20           induction       -0.2049576

HCPD

gameffects_neurosynth_hcpd <- merge(gameffects_FA_glasser_hcpd, neurosynth.terms, by = "label")

devpattern.neurosynth.hcpd <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_hcpd, dev.metric = "smooth.increase.offset", term = t)}) %>%
  arrange(desc(term.correlation))
devpattern.neurosynth.hcpd <- rbind(slice_max(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10))

devpattern.neurosynth.hcpd
##                      term term.correlation
## 1              expectancy        0.3933764
## 2  sentence_comprehension        0.3887129
## 3                 thought        0.3767177
## 4       cognitive_control        0.3477497
## 5     response_inhibition        0.3345571
## 6              monitoring        0.3258089
## 7                  recall        0.3033990
## 8      emotion_regulation        0.2986220
## 9                language        0.2927547
## 10              reasoning        0.2664190
## 11           multisensory       -0.3268444
## 12      visual_perception       -0.3132163
## 13     object_recognition       -0.3088589
## 14              induction       -0.2857823
## 15           localization       -0.2604849
## 16         mental_imagery       -0.2576200
## 17             adaptation       -0.2537276
## 18                   gaze       -0.2474086
## 19               fixation       -0.2473204
## 20             navigation       -0.2404549

Overlapping terms across datasets

merge(devpattern.neurosynth.pnc, devpattern.neurosynth.hcpd, by = "term") %>% arrange(desc(term.correlation.x)) %>% select(term)
##                   term
## 1           expectancy
## 2           monitoring
## 3    cognitive_control
## 4            reasoning
## 5  response_inhibition
## 6              thought
## 7               recall
## 8            induction
## 9                 gaze
## 10   visual_perception
## 11  object_recognition
ggplot(devpattern.neurosynth.pnc, aes(x = term.correlation, y = term, color = term.correlation)) +
  geom_segment(aes(y = reorder(term, term.correlation), yend = term, x = 0, xend =  
                     term.correlation), color = "#989898", linewidth = 0.2) +
  geom_point(size = 2.2, alpha = 0.85) +
  theme_light() +
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.29, 0.35), oob=squish) +
  labs(x = "\nCognitive Term-Thalamocortical Development Correlation", y = "NeuroSynth Term\n") +
  scale_x_continuous(limits = c(-0.3,0.42), breaks = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4)) +
  scale_y_discrete(
  labels = c("expectancy"=expression(bold(expectancy)),
  "monitoring"=expression(bold(monitoring)),
  "cognitive_control"=expression(bold(cognitive_control)),
  "reasoning"=expression(bold(reasoning)),
  "response_inhibition"=expression(bold(response_inhibition)),
  "thought"=expression(bold(thought)),
  "recall"=expression(bold(recall)),
  "induction"=expression(bold(induction)),
  "object_recognition"=expression(bold(object_recognition)),
  "visual_perception"=expression(bold(visual_perception)),
  "gaze"=expression(bold(gaze)),
  "induction"=expression(bold(induction), parse=TRUE))) +
   theme(
   legend.position = "none", 
   panel.grid.major.y = element_blank(),
   panel.grid.major.x = element_line(color = c("gray90")),
   axis.ticks.y = element_blank(),
   panel.border = element_blank(),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Neurosynth_Maturationdecoding.pdf", device = "pdf", dpi = 500, width = 3.35, height = 3.07)

Term overlap null analysis

#Create 10,000 spins of the HCPD age of maturation map for a neurosynth term overlap null analysis. Each spatially rotated map will be correlated with neurosynth term maps as above and the 10 most positively and negatively correlated terms will be identified. The number of overlapping terms with the empirical PNC terms will then be computed, creating a term overlap null based on spun maps for computing a pspin value.
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd) #complete data for spinning
x <- spin.data$smooth.increase.offset #specific cortical map to spin

#spin the empirical data 
perm.id <- perm.id.full
nroi = dim(perm.id)[1]  #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = array(NA,dim=c(nroi,nperm)) #initialize
for (r in 1:nperm) { #spin it
  for (i in 1:nroi) {
    x.perm[i,r] = x[perm.id[i,r]]
  }
}
x.spatialperm <- as.data.frame(x.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(x.perm))))
x.spatialperm$label <- spin.data$label

#compute neurosynth correlations for all spins and identify overlapping terms
compute.termoverlap <- function(perm){
  this.perm <- sprintf("perm%s", perm)
  nullmap.data <- x.spatialperm %>% select(label, all_of(this.perm))
  nullmap.neurosynth <- left_join(nullmap.data, neurosynth.terms, by = "label")
  nullmap.neurosynth <- left_join(nullmap.neurosynth, glasser.regions, by = "label")
  
  #null neurosynth correlations for this spin
  devpattern.neurosynth.null <- ldply(neurosynth.termlist, function(t){
    developmentmap.neurosynth.decoding(df = nullmap.neurosynth, dev.metric = this.perm, term = t)}) %>% 
    arrange(desc(term.correlation))
  
  #identify top and bottom most correlated terms for this spin
  devpattern.neurosynth.null <- rbind(slice_max(devpattern.neurosynth.null, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.null, order_by = term.correlation, n = 10))
  
  #compute overlap
  pnc.terms.empirical <- devpattern.neurosynth.pnc %>% select(term)
  null.terms <- devpattern.neurosynth.null %>% select(term)
  num.overlap <- length(intersect(pnc.terms$term, null.terms$term))
  
  return(num.overlap)
}

term.overlap.null <- lapply(1:nperm, compute.termoverlap)
term.overlap.null <- do.call(rbind, term.overlap.null) %>% as.data.frame %>% set_names("overlapping.terms")

write_rds(term.overlap.null, "/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/termoverlap_null.rds")
term.overlap.null <- readRDS("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/termoverlap_null.rds")

sprintf("In %s out of 10,000 spins, the term overlap number was equal to or greater than 11", sum(term.overlap.null$overlapping.terms >= 11))
## [1] "In 36 out of 10,000 spins, the term overlap number was equal to or greater than 11"
#calculate spin based p
sum(term.overlap.null$overlapping.terms >= 11) / 10000
## [1] 0.0036

Thalamocortical Connectivity Maturation is Organized by the Sensorimotor-Association Axis

gameffects_FA_glasser_pnc <- merge(gameffects_FA_glasser_pnc, (S.A.axis %>% select(tract, SA.axis)), by = "tract")
gameffects_FA_glasser_hcpd <- merge(gameffects_FA_glasser_hcpd, (S.A.axis %>% select(tract, SA.axis)), by = "tract")

Age-resolved developmental change plots

PNC

derivatives_FA_glasser_pnc$significant.derivative[derivatives_FA_glasser_pnc$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_pnc, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(size = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Derivativeplot_pnc.pdf", device = "pdf", dpi = 500, width = 2.3, height = 3)

HCPD

derivatives_FA_glasser_hcpd$significant.derivative[derivatives_FA_glasser_hcpd$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_hcpd, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(size = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

Thalamocortical connections mature at progressively older ages across the S-A axis

PNC

Spearman’s rho

cor.test(gameffects_FA_glasser_pnc$SA.axis, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_pnc$SA.axis and gameffects_FA_glasser_pnc$smooth.increase.offset
## S = 683025, p-value = 2.394e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4877183

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") #fdr-correct this for axis-based analyses below
## [1] 0.00075

Correlation plot

ggplot(gameffects_FA_glasser_pnc, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_pnc.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)

Brain map

Thalamocortical Age of Maturation

agemat.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
agemat.df <- left_join((spin.data %>% select(label)), agemat.df, by = "label")
agemat.df$cortex <- "cortex"
agemat.df <- rbind(agemat.df, data.frame(label = "rh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey
agemat.df <- rbind(agemat.df, data.frame(label = "lh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey

ggplot() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = smooth.increase.offset, colour=I("black"), size=I(.05))) +  
  theme_void()  + 
  scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint =  17.5, na.value = "gray93", limits = c(16, 19), oob = squish) + 
  new_scale_fill() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c(alpha("white", 0), "white")) +
  theme(legend.position = "none")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Ageofmaturation_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 5, height = 6.5) 

Sensorimotor-Association Axis

S.A.axis.plot <- S.A.axis[S.A.axis$tract %in% tractlist.pnc$tract,] %>% select(label, SA.axis)

ggplot() + 
  geom_brain(data = S.A.axis.plot, atlas = glasser, mapping = aes(fill = SA.axis, colour=I("black"), size=I(.05))) +  
  theme_void()  + 
  scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint =  180, na.value = "gray93") + 
   new_scale_fill() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c(alpha("white", 0), "white")) +
  theme(legend.position = "none")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 4.7, height = 6.2) 

HCPD

Spearman’s rho

cor.test(gameffects_FA_glasser_hcpd$SA.axis, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_hcpd$SA.axis and gameffects_FA_glasser_hcpd$smooth.increase.offset
## S = 367590, p-value = 2.934e-12
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5090025

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") #fdr-correct this for axis-based analyses below
## [1] 0.00155

Correlation plot

ggplot(gameffects_FA_glasser_hcpd, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ece4f2", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(10.7, 22), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_hcpd.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)

Thalamocortical maturation along anatomical axes

#Glasser parcel x,y,z coordinates for defining A-P, M-L, and D-V axes
load(file = "/cbica/projects/thalamocortical_development/Maps/parcellations/surface/hcp_mmp1.0.rda") #RDA file from the brainGraph package https://github.com/cwatson/brainGraph/blob/master/data/hcp_mmp1.0.rda with glasser region x, y, and z MNI coordinates. regions are in lh --> rh order
hcp_mmp1.0$x.mni[1:180]  <- hcp_mmp1.0$x.mni[1:180]*-1 #convert R-> L hemisphere coords into medial lateral coords
hcp_mmp1.0$label <- NA
hcp_mmp1.0$label[1:180] <- glasser.regions$label[181:360]
hcp_mmp1.0$label[181:360] <- glasser.regions$label[1:180]
hcp_mmp1.0 <- hcp_mmp1.0 %>% select(-hemi)
#C-Mt files for defining the group-level core-matrix gradient
CPt.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractaverage_CPt_primary.csv") %>% select("tract", "label", "orig_parcelname", "mean_CPt")
CPt.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractaverage_CPt_primary.csv")  %>% select("tract", "label", "orig_parcelname", "mean_CPt")
CPt.gradient <- merge(CPt.glasser.pnc, CPt.glasser.hcpd, by = c("label", "tract", "orig_parcelname"), all.x = TRUE, all.y = TRUE, suffixes = c("_pnc", "_hcpd"))
CPt.gradient$CPt <- CPt.gradient %>% select(mean_CPt_pnc, mean_CPt_hcpd) %>% rowMeans(na.rm = T)
CPt.gradient <- CPt.gradient %>% select(CPt, label)

hcp_mmp1.0 <- left_join(hcp_mmp1.0, CPt.gradient, by = "label")
#Function to compute the correlation between a development map and an anatomical axis (x,y,z) as well as test the significance of the difference between the developmental correlation with the anatomical axis versus the S-A axis
compute_axis_correlation <- function(df.dev, dev.metric, axis){
  df.dev.axes <- merge(df.dev, hcp_mmp1.0, by = "label")
  df.dev.axes.spin <- left_join(spin.parcels, df.dev.axes,  by = c("label", "tract", "orig_parcelname"))
  
  #S-A axis - development correlation (as computed above) for comparison of two overlapping correlations based on dependent groups
  SA.axis.cor <-  cor(df.dev.axes$SA.axis, df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman"))
  
  #Anatomical axis - development correlation
  anatomical.axis.cor <- cor(df.dev.axes[,axis], df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between anatomical axis and development metric
  anatomical.axis.pvalue <- perm.sphere.p(x = df.dev.axes.spin[,axis], y = df.dev.axes.spin[,dev.metric], perm.id = perm.id.full, corr.type = "spearman") #spin based p-value for the correlation
  
  #S-A axis - anatomical axis correlation
  SA.anatomical.cor <- cor(df.dev.axes$SA.axis, df.dev.axes[,axis], use = "complete.obs", method = c("spearman"))
  
  #Test for significant difference between two correlations based on dependent groups (i.e., correlations with an overlapping input)
  cocor.result <- cocor.pvalue <- cocor.dep.groups.overlap(
  r.jk = SA.axis.cor, 
  r.jh = anatomical.axis.cor, 
  r.kh = SA.anatomical.cor, 
  n = nrow(df.dev.axes), alternative = "two.sided", test = "hittner2003")
  
  cocor.pvalue <- cocor.result@hittner2003[["p.value"]]
  
  comparison.results <- data.frame(axis = axis, axis.cor = anatomical.axis.cor, axis.cor.pvalue = anatomical.axis.pvalue, SA.cor = SA.axis.cor, cocor.pvalue = cocor.pvalue)
  
  return(comparison.results)
  }

PNC

Anterior-posterior axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "y.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 y.mni 0.3101209          0.0621 0.4877183 0.0001552829

Medial-lateral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "x.mni")
##    axis   axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 x.mni 0.03024772         0.43765 0.4877183 3.832011e-08

Dorsal-ventral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "z.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 z.mni 0.1173292         0.26005 0.4877183 0.0002034701

Core-matrix gradient

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "CPt")
##   axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1  CPt 0.2897652          0.0287 0.4877183 0.0006272672

FDR-correct axis correlation ps

axis.cor.ps <- c(0.00075, 0.0621, 0.43765, 0.26005, 0.0287) #S-A, A-P, M-L, D-V, C-M
p.adjust(axis.cor.ps, method = c("fdr"))
## [1] 0.0037500 0.1035000 0.4376500 0.3250625 0.0717500

FDR-correct correlation comparison ps

cocor.ps <- c(0.0001552829, 3.832011e-08, 0.0002034701, 0.0006272672)
p.adjust(cocor.ps, method = c("fdr"))
## [1] 2.712935e-04 1.532804e-07 2.712935e-04 6.272672e-04
#Put above results into df for plotting
axis.results <- data.frame(Axis = factor(), corr = double())
axis.results <- axis.results %>% add_row(Axis = "S-A", corr = 0.4877183)
axis.results <- axis.results %>% add_row(Axis = "A-P", corr = 0.3101209)
axis.results <- axis.results %>% add_row(Axis = "C-M", corr = .2897652)
axis.results <- axis.results %>% add_row(Axis = "D-V", corr = 0.1173292)
axis.results <- axis.results %>% add_row(Axis = "M-L", corr = 0.03024772)
axis.results$Axis <- factor(axis.results$Axis, ordered = TRUE, levels = c("S-A", "A-P", "C-M", "D-V", "M-L"))

ggplot(axis.results, aes(x = Axis, y = corr)) + 
  geom_segment(aes(x = Axis, xend = Axis,  yend =  
                     corr, y = 0), linewidth = 0.2) +
  geom_point(aes(fill = Axis), size = 3.6, shape = 22, color = alpha("white", 0)) +
  scale_y_continuous(limits = c(-0.01, .5)) +
  labs(x="Neuroaxis") +
  labs(y="\nAxis Alignment") +
  theme_classic() +
  scale_fill_manual(values=c("#672975", "#9898BF", "#9898BF", "#9898BF", "#9898BF")) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.text.x = element_text(vjust = 0.7, angle = 20),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Anatomicalaxes_agematuration_comparison_lollyplot.pdf", device = "pdf", dpi = 500, width = 1.36, height = 1.05)

HCPD

Anterior-posterior axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "y.mni")
##    axis axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 y.mni 0.441284         0.02195 0.5090025     0.124339

Medial-lateral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "x.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 x.mni 0.1124713          0.3695 0.5090025 5.601018e-07

Dorsal-ventral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "z.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 z.mni 0.0814776         0.38565 0.5090025 1.390815e-05

Core-matrix gradient

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "CPt")
##   axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1  CPt 0.3809388          0.0163 0.5090025   0.02169413

Correct axis correlation ps

axis.cor.ps <- c(0.00155, 0.02195, 0.3695, 0.38565, 0.0163) #S-A, A-P, M-L, D-V, C-M
p.adjust(axis.cor.ps, method = c("fdr"))
## [1] 0.00775000 0.03658333 0.38565000 0.38565000 0.03658333

Correct correlation comparison ps

cocor.ps <- c(0.124339, 5.601018e-07, 1.390815e-05, 0.02169413)
p.adjust(cocor.ps, method = c("fdr"))
## [1] 1.243390e-01 2.240407e-06 2.781630e-05 2.892551e-02

Thalamocortical Connectivity Maturation Mirrors Timescales of Cortical Plasticity

Plasticity marker development data

#Fluctuation amplitude age of decrease onset data from Sydnor et al., 2023, Nat Neurosci https://www.nature.com/articles/s41593-023-01282-y 
boldamplitude.agedecrease.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/spatiotemp_dev_plasticity/cortical_maps/AgeofDeclineOnset_FirstNegDeriv.pscalar.nii") 

boldamplitude.dev <- data.frame(orig_parcelname = names(boldamplitude.agedecrease.cifti$Parcel), amplitude.agedecline = boldamplitude.agedecrease.cifti$data)
#T1/T2 ratio rate of developmental change from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
myelin.maxdev.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin.roc.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_mean_posterior_annualized_roc_myelin.pscalar.nii") #annualized rate of change

myelin.dev <- data.frame(orig_parcelname = names(myelin.maxdev.cifti$Parcel), myelin.agemaxdev = myelin.maxdev.cifti$data, myelin.annualroc = myelin.roc.cifti$data)
#Generate the E:I ratio-age slope map in Glasser HCP-MMP parcels from Zhang, Larsen et al., bioRxiv, https://www.biorxiv.org/content/10.1101/2023.06.22.546023v1.full

library(ciftiTools)

#Calculate E:I age slopes in Yan100 parcels
EI.group.ages <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/age_Yan.csv", header = F) %>% as.data.frame %>% set_names("age") %>% mutate(age = age/12) #average age of 29 age groups
EI.group.ratio <- t(read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_Yan.csv", header = F)) %>% as.data.frame() #average E:I ratio in all of the Yan100 parcels in the 29 age groups
Yan100.EI.ageslopes <- sapply(EI.group.ratio, function(x) lm(x ~ EI.group.ages$age)$coefficients[2]) %>% as.data.frame() %>% set_names("slope") #linear model coefficients for the association between age and E:I ratio in all Yan100 parcels

#Reorder Yan100 parcels to match the released version of the atlas
Yan100.region.reordering <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/Yanatlas_regionmapping.csv") #mapping of regions between version 
Yan100.EI.ageslopes$mapping <- Yan100.region.reordering$parcelnumber.released 
Yan100.EI.ageslopes <- Yan100.EI.ageslopes %>% arrange(mapping)

#Map parcel-level Yan100 data to fslr 32k vertex-level data
Yan100.densecifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/100Parcels_Yeo2011_17Networks.dscalar.nii")
Yan100.lh.vertices <- Yan100.densecifti$data$cortex_left %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.lh.vertices.slope <- left_join(Yan100.lh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
Yan100.rh.vertices <- Yan100.densecifti$data$cortex_right %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.rh.vertices.slope <- left_join(Yan100.rh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
lh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.leftcortex.csv", header = F)
rh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.rightcortex.csv", header = F)
Yan100.EI.ageslopes.densecifti <- as_cifti(cortexL = Yan100.lh.vertices.slope$slope, cortexR = Yan100.rh.vertices.slope$slope, cortexL_mwall = lh.medialwall$V1, cortexR_mwall = rh.medialwall$V1) #ciftify me captain
write_cifti(Yan100.EI.ageslopes.densecifti, "/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii")

#Parcellate the fslr 32k vertex-level E:I-age slope with the glasser HCP-MMP atlas
command = "-cifti-parcellate /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii /cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii -only-numeric"
ciftiTools::run_wb_cmd(command, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)

detach("package:ciftiTools", unload=TRUE)
EI.ageslope.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii") 
EIratio.dev <- data.frame(orig_parcelname = names(EI.ageslope.cifti$Parcel), EI.ageslope = EI.ageslope.cifti$data)
df.list <- list(boldamplitude.dev, myelin.dev, EIratio.dev) #dfs to merge
plasticity.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
plasticity.dev <- left_join(spin.parcels, plasticity.dev, by = "orig_parcelname")

Fluctuation amplitude decrease onset map

plasticity.dev$agedecline.protracted <- plasticity.dev$amplitude.agedecline
plasticity.dev$agedecline.protracted[plasticity.dev$agedecline.protracted == "NaN"] <- 23 #lets plot regions with most protracted functional expression of plasticity (amplitude never declines!)
snr.exclude <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/SNRmask_glasser360.csv") %>% filter(SNR.mask == 0) #but not those we excluded for low SNR
plasticity.dev$agedecline.protracted[plasticity.dev$label %in% snr.exclude$label] <- NA

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = agedecline.protracted, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c(alpha("#323280", 0.2), alpha("#323280", 0.4), alpha("#323280", 0.6), alpha("#323280", 0.9), alpha("#323280", 1)), limits = c(10, 17), oob = squish, na.value = "white") 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Boldamplitude_agedecline_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.7, height = 5.4) 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

T1/T2 annualized rate of increase map

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = myelin.annualroc, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c(alpha("#323280", 1), alpha("#323280", 0.8), alpha("#323280", 0.6), alpha("#323280", 0.4), alpha("#323280", 0.2)), limits = c(0.012, 0.016), oob = squish, na.value = "white") 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2ratio_annualizedroc_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.4, height = 5.13) 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

E:I ratio magnitude of developmental decline map

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = EI.ageslope, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c(alpha("#323280", 0.2), alpha("#323280", 0.4), alpha("#323280", 0.6), alpha("#323280", 0.9), alpha("#323280", 1)), limits = c(-0.04, -0.035), oob = squish, na.value = "white") 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIratio_decline_corticalmap.pdf", device = "pdf", dpi = 500, width = 5, height = 5.08) 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 6 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>,
## #   agedecline.protracted <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope agedecline.protracted
## 396         16.37678      0.008886078 -0.03484702              17.18593

Thalamocortical maturational patterns align to the spatiotemporal unfolding of plasticity readouts

PNC

plasticity.dev.pnc <- left_join(plasticity.dev, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))

Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline

Spearman’s rho

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$EI.ageslope, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$EI.ageslope
## S = 733425, p-value = 2.329e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4499176

Spatial permutation (spin) based p-value

EI.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.pnc, aes(x = EI.ageslope, y = smooth.increase.offset)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - myelin annualized rate of growth

Spearman’s rho

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$myelin.annualroc, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$myelin.annualroc
## S = 1929655, p-value = 3.141e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4472777

Spatial permutation (spin) based p-value

myelin.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.pnc, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset

Spearman’s rho

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$amplitude.agedecline, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$amplitude.agedecline
## S = 712074, p-value = 3.088e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3028346

Spatial permutation (spin) based p-value

amplitude.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman")

Correlation plot

ggplot(plasticity.dev.pnc, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(10, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

FDR-corrected p-values across correlations

plasticity.ps <- c(EI.thalamus.pnc, myelin.thalamus.pnc, amplitude.thalamus.pnc)
plasticity.ps
## [1] 0.00020 0.00090 0.02965
p.adjust(plasticity.ps, method = c("fdr"))
## [1] 0.00060 0.00135 0.02965

HCPD

plasticity.dev.hcpd <- left_join(plasticity.dev, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))

Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline

Spearman’s rho

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$EI.ageslope, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$EI.ageslope
## S = 409170, p-value = 9.57e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4534636

Spatial permutation (spin) based p-value

EI.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hcpd, aes(x = EI.ageslope, y = smooth.increase.offset, color = EI.ageslope)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - myelin annualized rate of growth

Spearman’s rho

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$myelin.annualroc, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$myelin.annualroc
## S = 1071660, p-value = 7.233e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4314374

Spatial permutation (spin) based p-value

myelin.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hcpd, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset

Spearman’s rho

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$amplitude.agedecline, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$amplitude.agedecline
## S = 378150, p-value = 7.236e-08
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4136816

Spatial permutation (spin) based p-value

amplitude.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hcpd, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(10, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

FDR-corrected p-values across correlations

plasticity.ps <- c(EI.thalamus.hcpd, myelin.thalamus.hcpd, amplitude.thalamus.hcpd)
plasticity.ps
## [1] 0.00010 0.00195 0.00890
p.adjust(plasticity.ps, method = c("fdr"))
## [1] 0.000300 0.002925 0.008900